/**
  MOMD project, Anyang Normal University, IMP-CAS
  \file ETFMath.h
  \class ETFMath
  \brief Math class, to provide some general math methods.
  \author SUN Yazhou, aisa.rabbit@163.com
  \date Created: 2020/07/09
  \date Last modified: 2020/10/05 by SUN Yazhou
  \copyright 2020 SUN Yazhou
  \copyright MOMD project, Anyang Normal University, IMP-CAS
*/

#ifndef ETFMath_h
#define ETFMath_h

#include <cmath>

class ETFMath{
public:
  ETFMath(){}
  virtual ~ETFMath(){}

  // physical functions
  static constexpr double Pi(){ return 3.14159265358979323846; }
  /// golden cut ratio
	static constexpr double Alpha(){ return 0.61803398874989484820458683436565; }
  /// rad per degree
	static constexpr double DEGREE(){ return 0.01745329251994329547; }
	static constexpr double Sqrt2(){ return 1.414213562373095048802; }
	static constexpr double Sqrt3(){ return 1.73205080756887719318; }
  static constexpr double uMeV(){ return 931.494061; } ///< u (C12=12.u) in MeV
  static constexpr double hbarc(){ return 197.32705; } ///< \retval hbar*c in MeV*fm
  static constexpr double FineStructureConstant(){ return 7.29735307964482e-3; } ///< 1./137.0359895
  static constexpr double Pion0Mass(){ return 134.9770; } ///< \retval pion0 mass in MeV
  static constexpr double PionPMMass(){ return 139.5706; } ///< \retval pion+- mass in MeV
  static constexpr double KaonPMMass(){ return 493.677; } ///< \retval pion+- mass in MeV
  static constexpr double c0(){ return 299.792458; } ///< \reteval light speed in mm/ns
  // CC = 0.321840605 = e0/(u0*c0*1E6) SI unit
  static constexpr double CC(){ return 0.32184060506023993; }

  // the z border of the effective magField (in mm)
	static const double kzMagIn, kzMagOut;
  static double zMagIn(){ return kzMagIn; }
  static double zMagOut(){ return kzMagOut; }

  /// \note r_global = (R.r_local)+r0
	static void rotate(const double *pIn, double *pOut, const double *angle);
  /// distance of a point to a 2-D straight line
  /// l: x = k*z+b or k*z-x+b=0
  static double DistanceToLine(double z, double x, double k, double b){
    return (k*z-x+b)/sqrt(1.+k*k);
  } // end member function DistanceToLine
  /// distance between 2 points in zOx plane
  static double Distance(double z0, double x0, double z1, double x1){
    return sqrt(sqr(z0-z1) + sqr(x0-x1));
  } // end member function DistanceToLine

  /// from Ek per u to velocity for a generic nucleon of mass uMeV
  static double EkPerUToBeta(double pEk){
    const double r = 1. + uMeV()/pEk;
    return sqrt(2.*r-1.)/r;
  }
  /// from Ek per u to momentum for a generic nucleon of mass uMeV
  static double EkPerUToP(double pEk){
    const double b = EkPerUToBeta(pEk);
    return b/sqrt(1.-b*b)*uMeV();
  }
  static double WoodsSaxon(double r, double a0, double r0){ return 1./(1.+exp((r-r0)/a0)); }

  /// sign of a number
  static double sign(double c){
    if(c >= 0.) return 1.;
    return -1.;
  }
  /// sign of a number
  static double sign(int c){ return sign(double(c)); }
  /// \retval (-)^m
  static double minus(int n){ return n & 1 ? -1 : 1; }
  static double sqr(double a){ return 0. == a ? a : a*a; }
  static double innerProd(const double *p0, const double *p1, int n = 3){
    double prod = 0.;
    for(int i = n; i--;) prod += p0[i]*p1[i];
    return prod;
  } // end member function innerProd
  static double norm(const double *p, int n = 3){ return sqrt(innerProd(p,p,n)); }
  /// \retval computes (a^2+b^2)^0.5 without desctructive underflow or overflow
  static double norm(double a, double b){
    a = fabs(a); b = fabs(b);
    if(a > b) return a*sqrt(1.+sqr(b/a));
    else return 0. == b ? 0. : b*sqrt(1.+sqr(a/b));
  } // end member function norm
  static double infNorm(int n, const double *f){
    double m = 0.;
    for(int i = n; i--;) if(fabs(f[i]) > m) m = fabs(f[i]);
    return m;
  } // end member function InfNorm
  // generic programming :)
  /// \retval: sum of square of each one in the parameter list
  static double Sum2(){ return 0.; }
  template <typename T, typename ...Args> // typename preferred, class deprecated
  static T Sum2(const T &v, const Args &...args){
    return v*v + Sum2(args...);
  }
  template <typename T>
  static bool Within(const T &v, const T &A, const T &B){
    if(v >= A && v <= B) return true;
    if(v >= B && v <= A) return true;
    return false;
  }
  static double ang(double k){
    return atan(k)/DEGREE();
  } // end member function
  static void peek(int *v, int len); ///< print array in an orderly fashion

  /// \retval calculate combination number
  /// returns the value ln[Gamma(xx)] for xx > 0.
  static double gammaln(double xx);
  // returns ln(n!)
  static double factln(int n);
  /// \retval n!
  static int Factorial(int x);
  /// \param n, m: n should not be smaller than m
  static int Binomial(int n, int m);
  /// beta function, Binomial coefficients with float arguments
  static double Beta(double z, double w);
  /// \retval n!!
  static int BiFactorial(int n);
  /// the incomplete gamma function gammap, gammpa(a,x)=P(a,x)=gamma(a,x)/Gamma(a)
  static double gammap(double a, double x);
  /// the incomplete gamma function gammap, gammpa(a,x)=P(a,x)=gamma(a,x)/Gamma(a)
  static double gammaq(double a, double x);
  /// returns the incomplete gamma function P(a,x) evaluated by its series
  /// representation. Optionally returns ln[Gamma(a)] as gln.
  static double gser(double a, double x, double *gln = nullptr);
  /// returns the incomplete gamma function P(a,x) evaluated by its continued fraction
  /// representation. Optionally returns ln[Gamma(a)] as gln.
  static double gcf(double a, double x, double *gln = nullptr);

  /// fit the track and assign slope and intercept
  /// \param np: number of points (z,x,r): DC hit position and the drift distance
  /// \return fitting chi2
  static double IterFit(const int np, const double *z, const double *x, const double *r,
    double &fK, double &fB);

  /// double rho(...): solve rho of track arc in the magnetic field
  /// \param (zo, xo): rotating center of the arc
  /// \param x2Arr: the two x2-s. x2 is the x coordinate of the exit trk from Mag
  /// the difference between x2-s from real trk and the calculated arc gives the pid estimate
  static double rho(double kin, double bin, double kout, double bout,
    double *x2Arr = nullptr, double *zo = nullptr, double *xo = nullptr);

	/// solve particle trajectory in uniform magnetic field, with only Mag boundry, exit track
	/// and target position known; returning the track radius of curvature in the magnetic field
	/// input unit: mm; output unit: mm
	/// x=ki z+bi, track after Target, before Mag; (zo, xo): rotating center of the arc
  /// \retval false: shit happens
	/// \return result: [0-7]: dtheta, rho, (ki, bi), (zo, xo), (x1, x2)
  // dtheta: deflection angle; rho, (zo, xo): defines the arc;
  // (ki, bi): defines the incident line; (x1, x2): x's of points of into and from magField
	static bool UniformMagneticSolution(double k1, double b1, double zMagOut,
		double zMagIn, double zTa, double xTa, double *result);

  static double BetaGamma(double beta);

  /// Given a set of data x[0...ndata-1] and y[0...ndata-1] with standard deviations
  /// sigy[0...ndata-1]. Fit them to a straight line by minimizing chi2. Returned
  /// are a and b, and their respective probable uncertainties.siga and sigb, the
  /// chi-square chi2 and the goodness-of-fit probability q (that the fit would
  /// have chi2 this large or larger under the model given by the fitting). If
  /// dy[0] < 0., then the standard deviations are assumed to be unavailable: q
  /// is returned as 1.0 and the normalization of chi2 is to unit standard deviation
  /// on all points.
  static void LinearFit(const double *x, const double *y, const double *dy, int ndata,
      double &a, double &b, double &siga, double &sigb, double &chi2, double &q);

  // ///////////////////// 3D TRACK PROJRECTION TRANSFORMATION ///////////////////////
	/// 3D track projections transformation, functions serving TAMWDCArray tracking.
	/// U+V->X tranformation: l: x=kz+b: slope
	static double kUV_X(double phi, double ku, double kv);
	/// U+V->X tranformation: l: x=kz+b: intercept
	static double bUV_X(double phi, double ku, double kv, double bu, double bv);
	/// U+V->Y tranformation: l: y=kz+b: slope
	static double kUV_Y(double phi, double ku, double kv);
	/// U+V->Y tranformation: l: y=kz+b: intercept
	static double bUV_Y(double phi, double ku, double kv, double bu, double bv);
	/// X+Y->U tranformation: l: xu=kzu+bu: slope
	static double kXY_U(double phi, double k1, double k2);
	/// X+Y->U tranformation: l: xu=kzu+bu: intercept
	static double bXY_U(double phi, double k1, double k2, double b1, double b2);
	/// X+Y->V tranformation: l: yv=kzv+bv: slope
	static double kXY_V(double phi, double k1, double k2);
	/// X+Y->V tranformation: l: yv=kzv+bv: intercept
	static double bXY_V(double phi, double k1, double k2, double b1, double b2);

  /// to calibrate absolute time to trigger time for PXI channels
  static double TimeToTrig(double t, unsigned bunchID){
    const double tt = t - (bunchID & 0x7FF) * 25.; // in ns
    return tt > 0. ? tt : tt + 51200.;
  } // end member function GetTimeToTrig
};

#endif
